package com.astrob.util;

import com.astrob.model.LonLat;
import com.nhc.model.CoordParmFileHandle;
import com.nhc.model.CoordSelfInfo;
import com.nhc.model.HCLocationInfo;
import com.nhc.model.LambertParmFileHandle;
import com.nhc.model.LambertSetInfo;
import com.nhc.model.PointAZInfo;
import com.nhc.model.TransParaSeven;
import com.nhc.model.TransParaSevenFileHandle;
import com.nhc.model.TwoPointsFileHandle;
import com.nhc.model.TwoPointsInfo;

import static java.lang.Math.pow;

/**
 * Created by colin on 2017/7/10.
 */


public class CoordChange {

    private static CoordChange instance = new CoordChange();

    public static CoordChange getInstance() {
        return instance;
    }

    private TransParaSeven transParaSeven = new TransParaSeven();


    public LonLat change2LL(String strLL, int coord_mode, int coord_format) {
        if (strLL == null) {
            return null;
        }

        int Zone = 0;//默认是0

        double pXcoord = 0;//lon
        double pYcoord = 0;//lat
        String UtmZone = "";//UtmZone + " " + strLon

        try {
            String[] names = strLL.split("\n");
            if (names == null || names.length != 2) {
                return null;
            }

            String[] strLat = names[0].split(" ");
            if (strLat == null || strLat.length != 2) {
                return null;
            }
            pYcoord = Double.parseDouble(strLat[0]);
            if (strLat[1].compareTo("S") == 0) {
                pYcoord = -pYcoord;
            }

            String[] strLon = names[1].split(" ");
            if (strLon == null) {
                return null;
            }

            if (strLon.length == 2) {
                pXcoord = Double.parseDouble(strLon[0]);
                if (strLon[1].compareTo("W") == 0) {
                    pXcoord = -pXcoord;
                }
            } else if (strLon.length == 3) {
                //"lat" -> "4426419.61 N"
                //"lon" -> "50S 442457.54 E"
                //"lon" -> "50SMK 442457.54 E"
                int len = strLon[0].length();
                if (len > 5) {
                    return null;
                } else if (len > 3) {
                    String strZone = strLon[0].substring(0, len - 3);
                    Zone = Integer.parseInt(strZone);
                    UtmZone = strLon[0].substring(len - 3);
                } else if (len >= 2) {
                    String strZone = strLon[0].substring(0, len - 1);
                    Zone = Integer.parseInt(strZone);
                    UtmZone = strLon[0].substring(len - 1);
                }

                pXcoord = Double.parseDouble(strLon[1]);
                if (strLon[2].compareTo("W") == 0) {
                    pXcoord = -pXcoord;
                }
            } else {
                return null;
            }
        } catch (Exception e) {
            return null;
        }

        //经纬度=>大地坐标=>经纬度 无法精确到小数点6位； CE版本大地坐标没有小数点
        return getLL(pXcoord, pYcoord, Zone, UtmZone, coord_mode, coord_format);
    }

    public LonLat change2LL(LonLat ll, int coord_mode, int coord_format) {
        if (ll == null) {
            return null;
        }

        //如果coord_format <= 2时

        return getLL(ll.GetLon(), ll.GetLat(), 0, "", coord_mode, coord_format);
    }

    //CFGUnitToLonLat--大地坐标=>经纬度   CFG2WGS
    public LonLat getLL(double pXcoord, double pYcoord, int Zone, String UtmZone, int coord_mode, int coord_format) {
        LonLat ll = new LonLat();

        transParaSeven = TransParaSevenFileHandle.getInstance().getListTransParaSeven(coord_mode);
        mgnDatumInfo_t CurDatum = new mgnDatumInfo_t();
        CurDatum.smaxis = transParaSeven.aAxis;//参考椭球长半轴
        CurDatum.flat = transParaSeven.flat;//参考椭球离心率
        CurDatum.delta[0] = -transParaSeven.DX;
        CurDatum.delta[1] = -transParaSeven.DY;
        CurDatum.delta[2] = -transParaSeven.DZ;
        CurDatum.scale = transParaSeven.scale;

        CurDatum.rot[0] = -transParaSeven.DXR;
        CurDatum.rot[1] = -transParaSeven.DYR;
        CurDatum.rot[2] = -transParaSeven.DZR;

        CenterLatitude = 0;//中央纬线
        CenterLongitude = 117;//中央经线
        EWOffset = 500000;//东西偏差
        NSOffset = 0;//南北偏差
        ProjectionScale = 1;//投影比例
        unitfactor = 1;
        
        double dMidLat = CenterLatitude * DEGTORAD;
        double dMidLon = CenterLongitude * DEGTORAD;

        mgnDGeodeticPos3D_t position = new mgnDGeodeticPos3D_t();
        position.Lon = pXcoord * DEGTORAD;
        position.Lat = pYcoord * DEGTORAD;

        XYcoord XYcoord = new XYcoord();
        XYcoord.pYcoord = pYcoord;
        XYcoord.pXcoord = pXcoord;
        UTMXY utmxy = new UTMXY();
        utmxy.pYcoord = pYcoord;
        utmxy.pXcoord = pXcoord;

        if (coord_format <= 2) {
            LoalToDatum(false, CurDatum, position);
            ll.SetLonLat(position.Lon, position.Lat);
        } else if (coord_format == 3) {
            //UTM
            if (UtmZone.compareTo("M") <= 0) {
                UtmZone = "S";
            } else {
                UtmZone = "N";
            }
            UTM2Geo(pXcoord, pYcoord, CurDatum, Zone, position, UtmZone);
            ll.SetLonLat(position.Lon, position.Lat);
        } else if (coord_format == 4) {
            //USER_TM
            CoordSelfInfo selfInfo = CoordParmFileHandle.getInstance().coordSelfInfo;
            CenterLatitude = selfInfo.ZLAT;//中央纬线
            CenterLongitude = selfInfo.ZLON;//中央经线
            EWOffset = selfInfo.EW;//东西偏差
            NSOffset = selfInfo.SN;//南北偏差
            ProjectionScale = selfInfo.TY;//投影比例

            dMidLat = CenterLatitude * DEGTORAD;
            dMidLon = CenterLongitude * DEGTORAD;
            UserXY2LLGeo2Tm(dMidLat, dMidLon, ProjectionScale,
                    position, XYcoord, EWOffset, NSOffset, unitfactor, CurDatum);
            ll.SetLonLat(position.Lon, position.Lat);
        } else if (coord_format == 5) {
            //MGRS
            int nGridScale = 1;
            int nOption = 3;
            MGRS2Geo(XYcoord, nGridScale, UtmZone, nOption, CurDatum, Zone, position);
            ll.SetLonLat(position.Lon, position.Lat);
        } else if (coord_format == 6) {
            //Lambert
            boolean OneStandPar = LambertParmFileHandle.getInstance().OneStandPar;
            LambertSetInfo lambert = null;
            if (OneStandPar) {
                lambert = LambertParmFileHandle.getInstance().lambert0;
            } else {
                lambert = LambertParmFileHandle.getInstance().lambert1;
            }
            double P1Lat = lambert.parallel1;//
            double P2Lat = lambert.parallel2;
            double lat1 = P1Lat * DEGTORAD;
            double lat2 = P2Lat * DEGTORAD;
            double orglat = lambert.ZLAT * DEGTORAD;
            double orglon = lambert.ZLON * DEGTORAD;
            UserXY2LLGeo2Lam(OneStandPar, orglat, orglon, lat1, lat2, lambert.TY,
                    position, XYcoord, lambert.EW, lambert.SN, lambert.utmc, CurDatum);
            ll.SetLonLat(position.Lon, position.Lat);
        } else if (coord_format == 7) {
            //Stereograph --立体投影
            int stereograph = LambertParmFileHandle.getInstance().stereograph;
            LambertSetInfo lambert = null;
            switch (stereograph) {
                default:
                case 0:
                    lambert = LambertParmFileHandle.getInstance().stereograph0;
                    break;
                case 1:
                    lambert = LambertParmFileHandle.getInstance().stereograph1;
                    break;
                case 2:
                    lambert = LambertParmFileHandle.getInstance().stereograph2;
                    break;
                case 3:
                    lambert = LambertParmFileHandle.getInstance().stereograph3;
                    break;

            }
            double orglat = lambert.ZLAT * DEGTORAD;
            double orglon = lambert.ZLON * DEGTORAD;
            UserXY2LLGeo2Ster(orglat, orglon, 10, lambert.TY,
                    position, XYcoord, lambert.EW, lambert.SN, lambert.utmc, CurDatum);
            ll.SetLonLat(position.Lon, position.Lat);
        } else if (coord_format == 8) {
            //Polyconic - 多圆锥投影
            LambertSetInfo lambert = LambertParmFileHandle.getInstance().polyconic;
            double orglat = lambert.ZLAT * DEGTORAD;
            double orglon = lambert.ZLON * DEGTORAD;
            UserXY2LLGeo2Poly(orglat, orglon, lambert.TY,
                    position, XYcoord, lambert.EW, lambert.SN, lambert.utmc, CurDatum);
            ll.SetLonLat(position.Lon, position.Lat);
        } else if (coord_format == 9) {
            //Oblique Mercator - 倾斜墨卡托投影
            double lat_o = 10;
            double lon_c = 117;
            double alpha_c = 200;
            double lat1 = 1;
            double lon1 = 1;
            double lat2 = 1;
            double lon2 = 1;

            int two_point = TwoPointsFileHandle.getInstance().mode;//
            if (two_point == 1) {
                TwoPointsInfo tpInfo = TwoPointsFileHandle.getInstance().mTwoPI;
                if (tpInfo != null) {
                    lat_o = tpInfo.ctnipt;
                    lat1 = tpInfo.lat_pt1;
                    lon1 = tpInfo.lon_pt1;
                    lat2 = tpInfo.lat_pt2;
                    lon2 = tpInfo.lon_pt2;
                    EWOffset = tpInfo.EW;//东西偏差
                    NSOffset = tpInfo.SN;//南北偏差
                    ProjectionScale = tpInfo.TY;//投影比例
                }
            } else {
                PointAZInfo pAZ = TwoPointsFileHandle.getInstance().mPAZ;
                if (pAZ != null) {
                    lat_o = pAZ.lat_o;
                    lon_c = pAZ.lon_c;
                    alpha_c = pAZ.alpha_c;
                    EWOffset = pAZ.EW;//东西偏差
                    NSOffset = pAZ.SN;//南北偏差
                    ProjectionScale = pAZ.TY;//投影比例
                }
            }
            lat_o = lat_o * DEGTORAD;
            lon_c = lon_c * DEGTORAD;
            alpha_c = alpha_c * DEGTORAD;
            lat1 = lat1 * DEGTORAD;
            lon1 = lon1 * DEGTORAD;
            lat2 = lat2 * DEGTORAD;
            lon2 = lon2 * DEGTORAD;
            UserXY2LLGeo2Obmer(two_point, lat_o, lon_c, alpha_c,
                    lat1, lon1, lat2, lon2, ProjectionScale,
                    position, XYcoord, EWOffset, NSOffset, unitfactor, CurDatum);
            ll.SetLonLat(position.Lon, position.Lat);
        } else if (coord_format >= 10 && coord_format <= 14) {
            //Philippines zone 1~5
            CoordSelfInfo selfInfo = CoordParmFileHandle.getInstance().getPhilippineszone(coord_format - 10);
            CenterLatitude = selfInfo.ZLAT;//中央纬线
            CenterLongitude = selfInfo.ZLON;//中央经线
            EWOffset = selfInfo.EW;//东西偏差
            NSOffset = selfInfo.SN;//南北偏差
            ProjectionScale = selfInfo.TY;//投影比例

            dMidLat = CenterLatitude * DEGTORAD;
            dMidLon = CenterLongitude * DEGTORAD;
            UserXY2LLGeo2Tm(dMidLat, dMidLon, ProjectionScale,
                    position, XYcoord, EWOffset, NSOffset, 1.0, CurDatum);
            ll.SetLonLat(position.Lon, position.Lat);
        } else if (coord_format == 15) {
            //Britn Grid
            EntBritCoord(position, CurDatum, utmxy);
            ll.SetLonLat(position.Lon, position.Lat);
        }

//        MGN_RAD2DEG
        double lat = MGN_RAD2DEG(ll.GetLat());
        double lon = MGN_RAD2DEG(ll.GetLon());
        ll.SetLonLat(lon, lat);
        return ll;
    }


    //--经纬度 => 大地坐标
    public HCLocationInfo getHCL(int mode, LonLat ll, int coord_format) {
        if (ll == null) {
            return null;
        }
        if (mode < 0 || mode >= SysUtils.COORD_0.length) {
            mode = 0;
        }

        if (mode == 0) {
            return new HCLocationInfo(SysUtils.COORD_0[mode], ll, "", coord_format);
        }
        
        CenterLatitude = 0;//中央纬线
        CenterLongitude = 117;//中央经线
        EWOffset = 500000;//东西偏差
        NSOffset = 0;//南北偏差
        ProjectionScale = 1;//投影比例
        unitfactor = 1;

        String UtmZone = "";
        double lon = ll.GetLon();
        double lat = ll.GetLat();
        transParaSeven = TransParaSevenFileHandle.getInstance().getListTransParaSeven(mode);

        mgnDatumInfo_t CurDatum = new mgnDatumInfo_t();
        CurDatum.smaxis = transParaSeven.aAxis;//参考椭球长半轴
        CurDatum.flat = transParaSeven.flat;//参考椭球离心率
        CurDatum.delta[0] = -transParaSeven.DX;
        CurDatum.delta[1] = -transParaSeven.DY;
        CurDatum.delta[2] = -transParaSeven.DZ;
        CurDatum.rot[0] = -transParaSeven.DXR;
        CurDatum.rot[1] = -transParaSeven.DYR;
        CurDatum.rot[2] = -transParaSeven.DZR;
        CurDatum.scale = transParaSeven.scale;

        mgnDGeodeticPos3D_t position = new mgnDGeodeticPos3D_t();
        position.Lat = lat * DEGTORAD;
        position.Lon = lon * DEGTORAD;
        UTMXY utmxy = new UTMXY();
        double dMidLat = CenterLatitude * DEGTORAD;
        double dMidLon = CenterLongitude * DEGTORAD;
        XYcoord XYcoord = new XYcoord();

        //mode--坐标系
        //coord_format --坐标投影
        switch (mode) {
            default: {
                if (coord_format <= 2) {
                    LoalToDatum(true, CurDatum, position);
                    ll.SetLonLat(position.Lon, position.Lat);
                    lon = MGN_RAD2DEG(position.Lon);
                    lat = MGN_RAD2DEG(position.Lat);
                } else if (coord_format == 3) {
                    //UTM
                    Geo2UTM(position, CurDatum, utmxy);
                    lon = utmxy.pXcoord;
                    lat = utmxy.pYcoord;
                    UtmZone = utmxy.pUtmZone + utmxy.UTMZoneStr;
                } else if (coord_format == 4) {
                    //USER_TM
                    CoordSelfInfo selfInfo = CoordParmFileHandle.getInstance().coordSelfInfo;
                    CenterLatitude = selfInfo.ZLAT;//中央纬线
                    CenterLongitude = selfInfo.ZLON;//中央经线
                    EWOffset = selfInfo.EW;//东西偏差
                    NSOffset = selfInfo.SN;//南北偏差
                    ProjectionScale = selfInfo.TY;//投影比例

                    dMidLat = CenterLatitude * DEGTORAD;
                    dMidLon = CenterLongitude * DEGTORAD;
                    UserLL2XYGeo2Tm(dMidLat, dMidLon, ProjectionScale,
                            position, XYcoord, EWOffset, NSOffset, unitfactor, CurDatum);
                    lon = XYcoord.pXcoord;
                    lat = XYcoord.pYcoord;
                } else if (coord_format == 5) {
                    //MGRS
                    int nGridScale = SysUtils.getInstance().mgrsaccuracy;//mgrsunit
                    int nOption = SysUtils.getInstance().mgrsgridtype;
                    Geo2MGRS(position, nGridScale, nOption, CurDatum, utmxy);
                    lon = utmxy.pXcoord;
                    lat = utmxy.pYcoord;
                    UtmZone = utmxy.pUtmZone + utmxy.UTMZoneStr;
                } else if (coord_format == 6) {
                    //Lambert
                    boolean OneStandPar = LambertParmFileHandle.getInstance().OneStandPar;
                    LambertSetInfo lambert = null;
                    if (OneStandPar) {
                        lambert = LambertParmFileHandle.getInstance().lambert0;
                    } else {
                        lambert = LambertParmFileHandle.getInstance().lambert1;
                    }
                    double lat1 = lambert.parallel1 * DEGTORAD;
                    double lat2 = lambert.parallel2 * DEGTORAD;
                    double orglat = lambert.ZLAT * DEGTORAD;
                    double orglon = lambert.ZLON * DEGTORAD;
                    UserLL2XYGeo2Lam(OneStandPar, orglat, orglon, lat1, lat2, lambert.TY,
                            position, XYcoord, lambert.EW, lambert.SN, lambert.utmc, CurDatum);
                    lon = XYcoord.pXcoord;
                    lat = XYcoord.pYcoord;

                } else if (coord_format == 7) {
                    //Stereograph --立体投影
                    int stereograph = LambertParmFileHandle.getInstance().stereograph;
                    LambertSetInfo lambert = null;
                    switch (stereograph) {
                        default:
                        case 0:
                            lambert = LambertParmFileHandle.getInstance().stereograph0;
                            break;
                        case 1:
                            lambert = LambertParmFileHandle.getInstance().stereograph1;
                            break;
                        case 2:
                            lambert = LambertParmFileHandle.getInstance().stereograph2;
                            break;
                        case 3:
                            lambert = LambertParmFileHandle.getInstance().stereograph3;
                            break;
                    }
                    double orglat = lambert.ZLAT * DEGTORAD;
                    double orglon = lambert.ZLON * DEGTORAD;
                    UserLL2XYGeo2Ster(orglat, orglon, 10, lambert.TY,
                            position, XYcoord, lambert.EW, lambert.SN, lambert.utmc, CurDatum);
                    lon = XYcoord.pXcoord;
                    lat = XYcoord.pYcoord;
                } else if (coord_format == 8) {
                    //Polyconic - 多圆锥投影
                    LambertSetInfo lambert = LambertParmFileHandle.getInstance().polyconic;
                    double orglat = lambert.ZLAT * DEGTORAD;
                    double orglon = lambert.ZLON * DEGTORAD;
                    UserLL2XYGeo2Poly(orglat, orglon, lambert.TY,
                            position, XYcoord, lambert.EW, lambert.SN, lambert.utmc, CurDatum);
                    lon = XYcoord.pXcoord;
                    lat = XYcoord.pYcoord;
                } else if (coord_format == 9) {
                    //Oblique Mercator - 倾斜墨卡托投影
                    double lat_o = 10;
                    double lon_c = 117;
                    double alpha_c = 200;
                    double lat1 = 1;
                    double lon1 = 1;
                    double lat2 = 1;
                    double lon2 = 1;

                    int two_point = TwoPointsFileHandle.getInstance().mode;//
                    if (two_point == 1) {
                        TwoPointsInfo tpInfo = TwoPointsFileHandle.getInstance().mTwoPI;
                        if (tpInfo != null) {
                            lat_o = tpInfo.ctnipt;
                            lat1 = tpInfo.lat_pt1;
                            lon1 = tpInfo.lon_pt1;
                            lat2 = tpInfo.lat_pt2;
                            lon2 = tpInfo.lon_pt2;
                            EWOffset = tpInfo.EW;//东西偏差
                            NSOffset = tpInfo.SN;//南北偏差
                            ProjectionScale = tpInfo.TY;//投影比例
                        }
                    } else {
                        PointAZInfo pAZ = TwoPointsFileHandle.getInstance().mPAZ;
                        if (pAZ != null) {
                            lat_o = pAZ.lat_o;
                            lon_c = pAZ.lon_c;
                            alpha_c = pAZ.alpha_c;
                            EWOffset = pAZ.EW;//东西偏差
                            NSOffset = pAZ.SN;//南北偏差
                            ProjectionScale = pAZ.TY;//投影比例
                        }
                    }

                    lat_o = lat_o * DEGTORAD;
                    lon_c = lon_c * DEGTORAD;
                    alpha_c = alpha_c * DEGTORAD;
                    lat1 = lat1 * DEGTORAD;
                    lon1 = lon1 * DEGTORAD;
                    lat2 = lat2 * DEGTORAD;
                    lon2 = lon2 * DEGTORAD;

                    UserLL2XYGeo2Obmer(two_point, lat_o, lon_c, alpha_c,
                            lat1, lon1, lat2, lon2, ProjectionScale,
                            position, XYcoord, EWOffset, NSOffset, unitfactor, CurDatum);
                    lon = XYcoord.pXcoord;
                    lat = XYcoord.pYcoord;
                } else if (coord_format >= 10 && coord_format <= 14) {
                    //Philippines zone 1~5
                    CoordSelfInfo selfInfo = CoordParmFileHandle.getInstance().getPhilippineszone(coord_format - 10);
                    CenterLatitude = selfInfo.ZLAT;//中央纬线
                    CenterLongitude = selfInfo.ZLON;//中央经线
                    EWOffset = selfInfo.EW;//东西偏差
                    NSOffset = selfInfo.SN;//南北偏差
                    ProjectionScale = selfInfo.TY;//投影比例

                    dMidLat = CenterLatitude * DEGTORAD;
                    dMidLon = CenterLongitude * DEGTORAD;
                    UserLL2XYGeo2Tm(dMidLat, dMidLon, ProjectionScale,
                            position, XYcoord, EWOffset, NSOffset, 1.0, CurDatum);
                    lon = XYcoord.pXcoord;
                    lat = XYcoord.pYcoord;

                } else if (coord_format == 15) {
                    //Britn Grid
                    GetBritCoord(position, CurDatum, utmxy);
                    lon = utmxy.pXcoord;
                    lat = utmxy.pYcoord;
                }
            }
            break;
            case 0:
                //WGS84
                if (coord_format > 2) {
                    coord_format = 2;
                }
            break;

        }
        return new HCLocationInfo(SysUtils.COORD_0[mode], new LonLat(lon, lat), UtmZone, coord_format);
    }


    private class mgnDGeodeticPos3D_t {
        public double Lat = 0;
        public double Lon = 0;
        public double Alt = 0;
    }

    private class mgnEcef_t {
        public double X = 0;
        public double Y = 0;
        public double Z = 0;
    }

    private class mgnDatumInfo_t {
        public String name = "";
        public double smaxis = 0;//参考椭球长半轴
        public double flat = 0;//参考椭球离心率
        public double[] delta = {0, 0, 0};//X,Y,Z
        public double[] rot = {0, 0, 0};
        public double scale = 0;
        public int flags = 7;
    }

    private class mgnDistArg_t {
        public double a_const = 0;
        public double b_const = 0;
        public double c_const = 0;
        public double d_const = 0;
        public double smaxis = 0;
        public double phi = 0;
    }

    private class XYcoord {
        public double pXcoord = 0;
        public double pYcoord = 0;
    }

    private class UTMXY {
        public double pXcoord = 0;
        public double pYcoord = 0;
        public int pUtmZone = 0;
        public String UTMZoneStr = "";
    }

    private class ENIndex {
        public int EastIndex = 0;
        public int NorthIndex = 0;
    }

    private final static double SEC2RAD = (1. / 206264.806247096);
    private final static double ESQ = 0.00669437999013;
    private final static double ERADUS = 6.378137e6;
    private final static double EPSQ = 0.00673949674227;
    private final static double SMI = 6356752.3142;
    private final static double TWOPI = (2.0 * Math.PI);
    private final static double UTM_PHI_0 = 0;

    //USER_TM----UTM
    private double CenterLatitude = 0;//中央纬线
    private double CenterLongitude = 117;//中央经线
    private int EWOffset = 500000;//东西偏差
    private int NSOffset = 0;//南北偏差
    private double ProjectionScale = 1;//投影比例
    private double unitfactor = 1;//通用米单位

    private final static double UTM_K_0 = 0.9996;
    private final static double RADTODEG = 180 / Math.PI;
    private final static double DEGTORAD = Math.PI / 180;
    private final static long UTM_FALSE_N = 10000000L;
    private final static long UTM_FALSE_E = 500000L;
    private final static int UTM_LAT_NUMBER = 20;

    private final static double HALFPI = Math.PI * 0.5;

    private void mgnLlh2Xyz(mgnDGeodeticPos3D_t Pos, mgnEcef_t XYZ) {
        double CosLat, SinLat; /**<  coslat, sinlat*/
        double Nuh, Nu;        /**<  nuh, nu;*/

        SinLat = Math.sin(Pos.Lat);
        CosLat = Math.cos(Pos.Lat);


        Nu = 1.0 - (ESQ * SinLat * SinLat);
        Nu = ERADUS / Math.sqrt(Nu);
        Nuh = Nu + (double) Pos.Alt;

        XYZ.X = Nuh * CosLat * Math.cos(Pos.Lon);
        XYZ.Y = Nuh * CosLat * Math.sin(Pos.Lon);
        XYZ.Z = (Nu * (1.0 - ESQ) + (double) Pos.Alt) * SinLat;
    }

    private void tnXyz2LlhLocal(mgnEcef_t pXyz, mgnDatumInfo_t Datum, mgnDGeodeticPos3D_t pPosition) {
        double R8Temp, pp, Costh, Sinth;
        double LocEradus, LocSmi, LocEsq, LocEpsq;

        LocEradus = Datum.smaxis; //loc_datum.smaxis;
        LocEsq = 2.0 * Datum.flat - Datum.flat * Datum.flat;
        LocEpsq = LocEsq / (1.0 - LocEsq);
        LocSmi = Datum.smaxis * (1.0 - Datum.flat);

        pp = Math.sqrt((pXyz.X) * (pXyz.X) + (pXyz.Y) * (pXyz.Y));
        R8Temp = Math.atan(((double) pXyz.Z * LocEradus) / (pp * LocSmi));
        Costh = Math.cos(R8Temp);
        Sinth = Math.sin(R8Temp);

        pPosition.Lat = Math.atan2(((double) pXyz.Z + LocEpsq * LocSmi * Sinth * Sinth * Sinth),
                (pp - LocEsq * LocEradus * Costh * Costh * Costh));

        pPosition.Lon = Math.atan2(pXyz.Y, pXyz.X);

        R8Temp = Math.sin(pPosition.Lat);
        R8Temp = 1.0 - LocEsq * R8Temp * R8Temp;
        pPosition.Alt = -LocEradus / Math.sqrt(R8Temp);
        Costh = Math.cos(pPosition.Lat);
        if (Costh != 0.0) {
            pPosition.Alt = pPosition.Alt + (pp / Costh);
        }
    }

    private void tnLlh2XyzLocal(mgnDGeodeticPos3D_t Position, mgnDatumInfo_t Datum, mgnEcef_t xyz) {
        double Coslat, Sinlat, nuh, nu;
        double LocEradus, LocEsq;

        LocEradus = Datum.smaxis;
        LocEsq = 2.0 * Datum.flat - Datum.flat * Datum.flat;

        Sinlat = Math.sin(Position.Lat);
        Coslat = Math.cos(Position.Lat);

        nu = 1.0 - (LocEsq * Sinlat * Sinlat);
        nu = LocEradus / Math.sqrt(nu);
        nuh = nu + (double) Position.Alt;

        xyz.X = nuh * Coslat * Math.cos(Position.Lon);
        xyz.Y = nuh * Coslat * Math.sin(Position.Lon);
        xyz.Z = (nu * (1.0 - LocEsq) + (double) Position.Alt) * Sinlat;
    }

    private void mgnXyz2Llh(mgnEcef_t XYZ, mgnDGeodeticPos3D_t Pos) {
        double r8temp, pp;   /**< double r8temp, pp */
        double Costh, Sinth; /**< double costh, sinth; */

        pp = Math.sqrt((XYZ.X) * (XYZ.X) + (XYZ.Y) * (XYZ.Y));
        double dTemp = Math.atan(((double) XYZ.Z * ERADUS) / (pp * SMI));

        r8temp = dTemp;
        Costh = Math.cos(r8temp);
        Sinth = Math.sin(r8temp);

        dTemp = Math.atan2(((double) XYZ.Z + EPSQ * SMI * Sinth * Sinth * Sinth),
                (pp - ESQ * ERADUS * Costh * Costh * Costh));

        Pos.Lat = dTemp;
        dTemp = Math.atan2(XYZ.Y, XYZ.X);
        Pos.Lon = dTemp;

        /**< HEIGHT IN METERS    */
        r8temp = Math.sin(Pos.Lat);
        r8temp = 1.0 - ESQ * r8temp * r8temp;
        Pos.Alt = -ERADUS / Math.sqrt(r8temp);
        Costh = Math.cos(Pos.Lat);

        if (Costh != 0.0) {
            Pos.Alt = Pos.Alt + (pp / Costh);
        }

    }

    private void LoalToDatum(boolean ToLocal, mgnDatumInfo_t CurDatum,
                             mgnDGeodeticPos3D_t position) {
        double rx, ry, rz;
        double sinx, siny, sinz, cosx, cosy, cosz;
        double[] r = {0, 0, 0, 0, 0, 0, 0, 0, 0};
        double scale = 1;
        mgnEcef_t xyz = new mgnEcef_t();
        mgnEcef_t xyz_tmp = new mgnEcef_t();

        mgnDGeodeticPos3D_t TmpPosition = position;
        //////////////////////////////////////////////////////////////////////////

		/* Convert rotation angles from seconds to radians */
        rx = CurDatum.rot[0] * SEC2RAD;
        ry = CurDatum.rot[1] * SEC2RAD;
        rz = CurDatum.rot[2] * SEC2RAD;

		/* Compute rotation matrix */
        sinx = Math.sin(rx);
        siny = Math.sin(ry);
        sinz = Math.sin(rz);
        cosx = Math.cos(rx);
        cosy = Math.cos(ry);
        cosz = Math.cos(rz);
        r[0] = cosy * cosz;
        r[1] = cosx * sinz + sinx * siny * cosz;
        r[2] = sinx * sinz - cosx * siny * cosz;
        r[3] = -cosy * sinz;
        r[4] = cosx * cosz - sinx * siny * sinz;
        r[5] = sinx * cosz + cosx * siny * sinz;
        r[6] = siny;
        r[7] = -sinx * cosy;
        r[8] = cosx * cosy;

        scale = 1. + CurDatum.scale * 1.e-6;

        if (ToLocal) {
            mgnLlh2Xyz(TmpPosition, xyz);

            xyz_tmp.X = r[0] * xyz.X + r[1] * xyz.Y + r[2] * xyz.Z;
            xyz_tmp.Y = r[3] * xyz.X + r[4] * xyz.Y + r[5] * xyz.Z;
            xyz_tmp.Z = r[6] * xyz.X + r[7] * xyz.Y + r[8] * xyz.Z;

		/* Scaling */
            xyz_tmp.X *= scale;
            xyz_tmp.Y *= scale;
            xyz_tmp.Z *= scale;

		/* Translation */
            xyz_tmp.X += CurDatum.delta[0];
            xyz_tmp.Y += CurDatum.delta[1];
            xyz_tmp.Z += CurDatum.delta[2];

            tnXyz2LlhLocal(xyz_tmp, CurDatum, TmpPosition);

        } else {
            tnLlh2XyzLocal(TmpPosition, CurDatum, xyz); //llh2xyz_local(*rlat, *rlon, alt, xyz, datum);

		/* Translation */
            xyz.X -= CurDatum.delta[0];
            xyz.Y -= CurDatum.delta[1];
            xyz.Z -= CurDatum.delta[2];

		/* Scaling */
            xyz.X /= scale;
            xyz.Y /= scale;
            xyz.Z /= scale;

		/* Rotation */
            xyz_tmp.X = r[0] * xyz.X + r[3] * xyz.Y + r[6] * xyz.Z;
            xyz_tmp.Y = r[1] * xyz.X + r[4] * xyz.Y + r[7] * xyz.Z;
            xyz_tmp.Z = r[2] * xyz.X + r[5] * xyz.Y + r[8] * xyz.Z;

            mgnXyz2Llh(xyz_tmp, TmpPosition); /** xyz2llh(xyz_tmp, rlat, rlon, &alt) */
        }

        position = TmpPosition;

    }

    private boolean UserLL2XYGeo2Tm(double lat, double lon, double k_0,
                                    mgnDGeodeticPos3D_t tempPos,
                                    XYcoord XYcoord,
                                    double false_e, double false_n, double unitfactor,
                                    mgnDatumInfo_t localDatum) {

        LoalToDatum(true, localDatum, tempPos);

        Geo2Tm(true, lat, lon, k_0, tempPos, XYcoord, localDatum);

        XYcoord.pXcoord = (XYcoord.pXcoord + false_e) / unitfactor;
        XYcoord.pYcoord = (XYcoord.pYcoord + false_n) / unitfactor;

        if ((XYcoord.pXcoord < -9999999.0) || (XYcoord.pXcoord > 99999999.0) ||
                (XYcoord.pYcoord < -9999999.0) || (XYcoord.pYcoord > 99999999.0)) {
            return false;
        }

        return true;
    }

    private void Geo2Tm(boolean Geo2Tm,
                        double lat, double lon, double k_0,
                        mgnDGeodeticPos3D_t tempPos,
                        XYcoord XYcoord,
                        mgnDatumInfo_t mgnDatum) {

        double mgnlat = lat;
        double mgnlon = lon;
        double mgnk_0 = k_0;

        mgnGeo2Tm(Geo2Tm, mgnlat, mgnlon, mgnk_0, tempPos, XYcoord, mgnDatum);
    }

    private void mgnGeo2Tm(boolean Geo2Tm, double lat, double lon, double k_0, mgnDGeodeticPos3D_t tempPos,
                           XYcoord XYcoord, mgnDatumInfo_t pDatumBlock) {

        double e_2, e_4, e_6;
        double one_minus_e2, e_prime_2;
        double e_const, f_const, m_0;
        double n_temp;
        double tan_phi, c_temp, a_temp;
        double a_temp_2, a_temp_3, a_temp_4, a_temp_5, a_temp_6;
        double mu, m_temp, e_1, phi__1;
        double cos_phi_1, c_1, tan_phi_1, n_1, temp_double;
        double r_1, d_temp, d_2, d_3, d_4, d_5, d_6, c_1_2, t_1_2;
        double e_1_2, e_1_3, e_1_4;
        double t_1, t_temp;
        double t_temp_2;
        double two_mu, cos_2_mu;
        double temp_1;
        mgnDistArg_t m_args = new mgnDistArg_t();

        e_2 = 1.0 - (1.0 - pDatumBlock.flat) * (1.0 - pDatumBlock.flat);
        one_minus_e2 = 1.0 - e_2;
        e_4 = e_2 * e_2;
        e_6 = e_4 * e_2;
        e_prime_2 = e_2 / one_minus_e2;
        m_args.a_const = -e_2 * 0.25;
        m_args.b_const = -e_4 * 0.046875;
        m_args.c_const = -e_6 * 0.01953125;
        m_args.d_const = 1.0 + m_args.a_const + m_args.b_const + m_args.c_const;
        m_args.smaxis = pDatumBlock.smaxis;

        m_args.phi = lat;
        m_0 = tnMDist(m_args);

        if (Geo2Tm) {

            temp_1 = Math.sin(tempPos.Lat);
            n_temp = pDatumBlock.smaxis / Math.sqrt(1.0 - e_2 * temp_1 * temp_1);
            tan_phi = Math.tan(tempPos.Lat);
            t_temp = tan_phi * tan_phi;
            t_temp_2 = t_temp * t_temp;
            temp_1 = Math.cos(tempPos.Lat);
            c_temp = e_prime_2 * temp_1 * temp_1;

            a_temp = tempPos.Lon - lon;
            if (a_temp > Math.PI) {
                a_temp = a_temp - TWOPI;
            } else if (a_temp < -Math.PI) {
                a_temp = a_temp + TWOPI;
            }
            a_temp = a_temp * temp_1;
            a_temp_2 = a_temp * a_temp;
            a_temp_3 = a_temp_2 * a_temp;
            a_temp_4 = a_temp_2 * a_temp_2;
            a_temp_5 = a_temp_4 * a_temp;
            a_temp_6 = a_temp_3 * a_temp_3;

            m_args.phi = tempPos.Lat;
            m_temp = tnMDist(m_args);

            XYcoord.pXcoord = k_0 * n_temp * (a_temp + (1.0 - t_temp + c_temp) * 0.166666667 * a_temp_3 +
                    (5.0 - 18.0 * t_temp + t_temp_2 + 72.0 * c_temp -
                            58.0 * e_prime_2) * 0.00833333333 * a_temp_5);

            XYcoord.pYcoord = k_0 * (m_temp - m_0 + n_temp * tan_phi * (0.5 * a_temp_2 +
                    (5.0 - t_temp + 9.0 * c_temp + 4.0 * c_temp * c_temp) * 0.0416666667 * a_temp_4 +
                    (61.0 - 58.0 * t_temp + t_temp_2 + 600.0 * c_temp - 330.0 * e_prime_2) *
                            0.00138888889 * a_temp_6));

        } else {
            mu = ((XYcoord.pYcoord / k_0) + m_0) / (pDatumBlock.smaxis * m_args.d_const);
            temp_double = Math.sqrt(one_minus_e2);
            e_1 = (1.0 - temp_double) / (1.0 + temp_double);
            e_1_2 = e_1 * e_1;
            e_1_3 = e_1_2 * e_1;
            e_1_4 = e_1_2 * e_1_2;

		/* footprint latitude, lat of y at central meridian
        */
            two_mu = 2.0 * mu;
            cos_2_mu = Math.cos(two_mu);
            e_const = e_1_3 * 2.416666666666667;
            f_const = e_1_4 * 12.0078125;

		/* footprint latitude, lat of y at central meridian
        */
            phi__1 = mu + Math.sin(two_mu) * (e_1 * 1.5 - e_const + cos_2_mu * (e_1_2 * 2.625
                    - f_const + cos_2_mu * (e_const * 2.603448276
                    + cos_2_mu * f_const * 1.427456083)));

            cos_phi_1 = Math.cos(phi__1);
            c_1 = e_prime_2 * cos_phi_1 * cos_phi_1;
            tan_phi_1 = Math.tan(phi__1);
            t_1 = tan_phi_1 * tan_phi_1;
            temp_double = Math.sin(phi__1);
            temp_double = Math.sqrt(1.0 - e_2 * temp_double * temp_double);
            n_1 = pDatumBlock.smaxis / temp_double;
            r_1 = pDatumBlock.smaxis * one_minus_e2 /
                    (temp_double * temp_double * temp_double);
            d_temp = XYcoord.pXcoord / (n_1 * k_0);
            d_2 = d_temp * d_temp;
            d_3 = d_2 * d_temp;
            d_4 = d_2 * d_2;
            d_5 = d_4 * d_temp;
            d_6 = d_3 * d_3;
            c_1_2 = c_1 * c_1;
            t_1_2 = t_1 * t_1;

            tempPos.Lat = phi__1 - (n_1 * tan_phi_1 / r_1) * (0.5 * d_2 - 0.0416666667 * d_4 *
                    (5.0 + 3.0 * t_1 + 10.0 * c_1 - 4.0 * c_1_2 - 9.0 * e_prime_2) +
                    (61.0 + 90.0 * t_1 + 298.0 * c_1 + 45.0 * t_1_2 - 252.0 * e_prime_2 -
                            3.0 * c_1_2) * 0.00138888889 * d_6);

            tempPos.Lon = lon + (d_temp - (1.0 + 2.0 * t_1 + c_1) * 0.16666667 * d_3 +
                    (5.0 - 2.0 * c_1 + 28.0 * t_1 - 3.0 * c_1_2 + 8.0 * e_prime_2 + 24.0 * t_1_2) *
                            0.00833333333 * d_5) / cos_phi_1;
        }
    }

    private double tnMDist(mgnDistArg_t args) {
        double two_phi, cos_2_phi, temp_2, temp_3;

        two_phi = 2.0 * args.phi;
        cos_2_phi = Math.cos(two_phi);
        temp_2 = 1.5 * args.a_const + 2.0 * args.b_const + 1.666666666666667 * args.c_const;
        temp_3 = 2.5 * args.b_const + 4.5 * args.c_const - 2.3333333333333333 * args.c_const * cos_2_phi;
        return (args.smaxis * (args.phi * args.d_const + Math.sin(two_phi) * (temp_2 - cos_2_phi * temp_3)));
    }


    //坐标显示格式UTM+自定义
    //Geo2UTM =>  UtmZone=50  UTMZoneStr=S  Xcoord=4426494.58N	Ycoord=442608.61E
    private boolean Geo2UTM(mgnDGeodeticPos3D_t tempPos, mgnDatumInfo_t localDatum, UTMXY utmxy) {
        LoalToDatum(true, localDatum, tempPos);
        return mgnGetUTMCoord(tempPos, localDatum, utmxy);
    }

    private boolean mgnGetUTMCoord(mgnDGeodeticPos3D_t Pos, mgnDatumInfo_t pDatumBlock, UTMXY utmxy) {
        /**< UTM valid between 80 degrees south and 84 degrees north plus 30 minutes overlap */
        if (Pos.Lat > 1.474803218 || Pos.Lat < -1.404990048)
            return false;

        int zone = mgnUTMZone(Pos.Lat, Pos.Lon);

        XYcoord XYcoord = new XYcoord();
        mgnGeo2Tm(true, 0, UTMcnmer(zone), UTM_K_0, Pos, XYcoord, pDatumBlock);

        if (XYcoord.pYcoord < 0.0)
            XYcoord.pYcoord += UTM_FALSE_N;

        utmxy.pXcoord = XYcoord.pXcoord + UTM_FALSE_E + 0.5;
        utmxy.pYcoord = XYcoord.pYcoord + 0.5;
        utmxy.pUtmZone = zone;

        int band = mgnUTMBand(Pos.Lat, Pos.Lon);
        utmxy.UTMZoneStr = mgnGetUTMLatBand(band);

        return true;

    }

    private int mgnUTMZone(double phi, double lam) {
        int zone;
        int zone_phi, zone_lam;

        lam = lam * RADTODEG;
        if (lam < 0.0)
            lam += 360.0;

        zone_phi = (int) (phi * RADTODEG);
        zone_lam = (int) lam;

        zone = (zone_lam / 6) + 31;
        if (zone > 60)
            zone -= 60;

        /**< Norwegian zone exceptions */
        if (zone_phi >= 56 && zone_phi < 64) {
            if (zone_lam >= 3 && zone_lam < 12) {
                zone = 32;
            }
        } else if (zone_phi >= 72) {
            if (zone_lam >= 0 && zone_lam < 9) {
                zone = 31;
            } else if (zone_lam >= 9 && zone_lam < 21) {
                zone = 33;
            } else if (zone_lam >= 21 && zone_lam < 33) {
                zone = 35;
            } else if (zone_lam >= 33 && zone_lam < 42) {
                zone = 37;
            }
        }
        return zone;
    }

    double UTMcnmer(int zone) {
        int temp_int = 0;

        if (zone > 30) {
            temp_int = ((zone - 31) * 6) + 3;
        } else {
            temp_int = -(((30 - zone) * 6) + 3);
        }
        return ((double) temp_int * DEGTORAD);
    }

    int mgnUTMBand(double phi, double lam) {
        phi = phi * RADTODEG;
        if (phi < -80)
            return -1;
        if (phi > 84)
            return -1;
        if (phi > 80 && phi <= 84)
            return UTM_LAT_NUMBER - 1;
        int band = (int) ((phi + 80) / 8);
        return band;
    }

    public static final String[] pUTMLatBand = {
            "C", "D", "E", "F", "G",
            "H", "J", "K", "L", "M",
            "N", "P", "Q", "R", "S",
            "T", "U", "V", "W", "X"};

    String mgnGetUTMLatBand(int Index) {
        if ((Index < 0) || (Index >= UTM_LAT_NUMBER)) {
            return "";
        }
        return pUTMLatBand[Index];
    }


    //其他坐标投影

    //mgnDGeodeticPos3D_t tempPos -- 输入的经纬度 double dLat, double dLon
    //UTMXY utmxy  -- 输出 pXcoord+pYcoord+pUtmZone+UTMZoneStr
    //--- int &pUtmZone, double &pXcoord, double &pYcoord, char *MGRSZoneStr
    //mgnDatumInfo_t localDatum -- 输入 椭球的7参数

    //MGRS
    boolean Geo2MGRS(mgnDGeodeticPos3D_t tempPos, int nGridScale, int nOption, mgnDatumInfo_t localDatum, UTMXY utmxy) {
        return mgnGetMgrsCoord(tempPos, nGridScale, nOption, localDatum, utmxy);
    }

    double MGN_DEG2RAD(double deg) {
        return deg * DEGTORAD;
    }

    double MGN_RAD2DEG(double rad) {
        return rad * RADTODEG;
    }

    private final static double HALF_MET_RAD = 78.5335e-9;

    public static final String MGRS_EAST_LETTER[][] = new String[][]{
            {"S", "T", "U", "V", "W", "X", "Y", "Z"},
            {"A", "B", "C", "D", "E", "F", "G", "H"},
            {"J", "K", "L", "M", "N", "P", "Q", "R"}
    };

    public static final String MGRS_EAST_LETTER2[] = {
            "S", "T", "U", "V", "W", "X", "Y", "Z",
            "A", "B", "C", "D", "E", "F", "G", "H",
            "J", "K", "L", "M", "N", "P", "Q", "R"
    };

    public static final String[] MGRS_GRID_LETTER = {
            "A", "B", "C", "D", "E",
            "F", "G", "H", "J", "K",
            "L", "M", "N", "P", "Q",
            "R", "S", "T", "U", "V",
            "W", "X", "Y", "Z"};

    boolean mgnGetMgrsCoord(mgnDGeodeticPos3D_t Pos, int MgrsRef, int mgrstopt, mgnDatumInfo_t pDatumBlock, UTMXY utmxy) {
        int TableOption;
        int Index, EastIndex, NorthIndex;
        int GridPower;
        long LongXcoord, LongYcoord;

        if (Pos.Lat > MGN_DEG2RAD(84.5) || Pos.Lat < MGN_DEG2RAD(-80.5))
            return false;

        int zone = mgnUTMZone(Pos.Lat, Pos.Lon);
        XYcoord XYcoord = new XYcoord();
        mgnGeo2Tm(true, 0, UTMcnmer(zone), UTM_K_0, Pos, XYcoord, pDatumBlock);

        if (Pos.Lat >= 0.0) {
            Index = ((int) ((Pos.Lat + HALF_MET_RAD) * RADTODEG) / 8) + 10;
            if (Index > 19)
                Index = 19;
        } else {
            Index = 9 - ((int) (-(Pos.Lat + HALF_MET_RAD) * RADTODEG) / 8);
            if (Index < 0)
                Index = 0;
        }

        utmxy.UTMZoneStr = mgnGetUTMLatBand(Index);
        EastIndex = (int) (((long) XYcoord.pXcoord + 400000L) / 100000L);
        utmxy.UTMZoneStr += MGRS_EAST_LETTER[zone % 3][EastIndex];
        if (XYcoord.pYcoord < 0.0)
            XYcoord.pYcoord += 1.0e+7;

        while (XYcoord.pYcoord >= 2.0e+6) {
            XYcoord.pYcoord -= 2.0e+6;
        }
        NorthIndex = (int) (XYcoord.pYcoord * 1.0e-5);

        switch (mgrstopt) {
            case 0: /**< MGN_MGRS_TP */
                TableOption = 1;//MGRS_TABLE_PRELIM;
                break;
            case 1: /**< MGRS_TE */
                TableOption = 2;//MGRS_TABLE_EDONE;
                break;
            default:
                TableOption = 4;//MGRS_TABLE_ALL;
        }


        int zz = zone % 2;

        if ((pDatumBlock.flags & TableOption) != 0) {
            if (zz == 0) NorthIndex += 5;
        } else {
            if (zz != 0) {
                NorthIndex += 10;
            } else {
                NorthIndex += 15;
            }
        }

        if (NorthIndex > 19)
            NorthIndex -= 20;
        utmxy.UTMZoneStr += MGRS_GRID_LETTER[NorthIndex];
        XYcoord.pXcoord += UTM_FALSE_E;

        while (XYcoord.pXcoord >= 1.0e+5) {
            XYcoord.pXcoord -= 1.0e+5;
        }

        while (XYcoord.pYcoord >= 1.0e+5) {
            XYcoord.pYcoord -= 1.0e+5;
        }

        GridPower = MgrsRef;

        LongXcoord = (long) ((XYcoord.pXcoord * pow(10.0, -GridPower)) + 0.5);
        LongYcoord = (long) ((XYcoord.pYcoord * pow(10.0, -GridPower)) + 0.5);

        if (LongXcoord == pow(10.0, 5 - GridPower))
            --LongXcoord;

        if (LongYcoord == pow(10.0, 5 - GridPower))
            --LongYcoord;

        utmxy.pUtmZone = zone;
        utmxy.pXcoord = (double) LongXcoord;
        utmxy.pYcoord = (double) LongYcoord;
        return true;
    }

    //Lambert
    private boolean UserLL2XYGeo2Lam(boolean OneStandPar, double lat, double lon, double lat1, double lat2,
                                     double k_0, mgnDGeodeticPos3D_t tempPos,
                                     XYcoord XYcoord, double false_e, double false_n, double unitfactor,
                                     mgnDatumInfo_t localDatum) {

        LoalToDatum(true, localDatum, tempPos);

        Geo2Lam(true, OneStandPar, lat, lon, lat1, lat2, k_0, tempPos, XYcoord, localDatum);

        XYcoord.pXcoord = (XYcoord.pXcoord + false_e) / unitfactor;
        XYcoord.pYcoord = (XYcoord.pYcoord + false_n) / unitfactor;

        if ((XYcoord.pXcoord < -9999999.0) || (XYcoord.pXcoord > 99999999.0) ||
                (XYcoord.pYcoord < -9999999.0) || (XYcoord.pYcoord > 99999999.0)) {
            return false;
        }


        return true;
    }

    private void Geo2Lam(boolean Geo2Lam, boolean OneStandPar,
                         double lat, double lon, double lat1, double lat2, double k_0,
                         mgnDGeodeticPos3D_t tempPos,
                         XYcoord XYcoord,
                         mgnDatumInfo_t mgnDatum) {

        double mgnPhi0 = lat;
        double mgnLamb0 = lon;
        double mgnPhi1 = lat1;
        double mgnPhi2 = lat2;
        double mgnk_0 = k_0;

        mgnGeo2Lam(Geo2Lam, OneStandPar, mgnPhi0, mgnLamb0, mgnPhi1, mgnPhi2, mgnk_0, tempPos, XYcoord, mgnDatum);
    }

    double tnTfunc(double phi, double ecc) {
        double SinPhi, eSinPhi;
        SinPhi = Math.sin(phi);
        eSinPhi = ecc * SinPhi;
        return (Math.sqrt(((1.0 - SinPhi) / (1.0 + SinPhi))
                * pow((1.0 + eSinPhi) / (1.0 - eSinPhi), ecc)));
    }

    double tnMfunc(double phi, double e_2) {
        double SinPhi;
        SinPhi = Math.sin(phi);
        return (Math.cos(phi) / Math.sqrt(1.0 - e_2 * SinPhi * SinPhi));
    }

    double LambDiff(double lamb, double lamb_0) {
        return (InitPIBounds(lamb - lamb_0));
    }

    double InitPIBounds(double angle) {
        if (angle > Math.PI) {
            return (angle - TWOPI);
        } else if (angle < -Math.PI) {
            return (angle + TWOPI);
        } else {
            return (angle);
        }
    }

    double Chi2Phi(double LocalChi, mgnDatumInfo_t pCurDatum) {

        double twoChi, cos2Chi;
        double aConst, bConst, cConst;
        double E_4, E_6, E_8;
        double TmpDouble, E_2;

        TmpDouble = 1.0 - pCurDatum.flat;
        E_2 = 1.0 - TmpDouble * TmpDouble;
        E_4 = E_2 * E_2;
        E_6 = E_4 * E_2;
        E_8 = E_4 * E_4;

        /**< conformal to geodetic lat */
        twoChi = 2.0 * LocalChi;
        cos2Chi = Math.cos(twoChi);
        aConst = E_4 * 0.2083333333333;
        bConst = E_6 * 0.025;
        cConst = E_8 * 0.03621031746;

        return (LocalChi + Math.sin(twoChi) * (0.5 * E_2 + aConst + bConst - cConst +
                cos2Chi * (aConst * 1.4 + bConst * 9.666666666666667 +
                        cConst * 0.9575342466 + cos2Chi * (bConst * 9.333333333333333 +
                        cConst * 7.989041096 + cos2Chi * cConst * 5.861643836))));
    }

    private void mgnGeo2Lam(boolean Geo2Lamb, boolean OneStandPar, double Phi0, double Lamb0,
                            double Phi1, double Phi2, double K0, mgnDGeodeticPos3D_t tempPos,
                            XYcoord XYcoord, mgnDatumInfo_t pCurDatum) {
        double TempDouble;
        double OneMinusE2, Ecc;
        double E_2, E_4, E_6;
        double Ntemp; //n_temp;
        double TempDouble1; //temp_double_1;
        double AScaled, T0, T1, T2, M1, M2, Ftemp, Rho0;
        double Theta, Ttemp, Rho, SignN, Chi;

        TempDouble = 1.0 - pCurDatum.flat;
        E_2 = 1.0 - TempDouble * TempDouble;
        OneMinusE2 = 1.0 - E_2;
        Ecc = Math.sqrt(E_2);
        E_4 = E_2 * E_2;
        E_6 = E_4 * E_2;

        AScaled = pCurDatum.smaxis * K0;

        T0 = tnTfunc(Phi0, Ecc);
        T1 = tnTfunc(Phi1, Ecc);
        M1 = tnMfunc(Phi1, E_2);

        if (OneStandPar) {
            Ntemp = Math.sin(Phi1);
        } else {
            T2 = tnTfunc(Phi2, Ecc);
            M2 = tnMfunc(Phi2, E_2);
            TempDouble1 = Math.log(M1) - Math.log(M2);
            TempDouble = Math.log(T1) - Math.log(T2);
            Ntemp = TempDouble1 / TempDouble;
            // equivalent to Richardus' sin(Phi0), i.e. sine of calculated standard parallel
        }

        Ftemp = M1 / (Ntemp * pow(T1, Ntemp));

        Rho0 = AScaled * Ftemp * pow(T0, Ntemp);

        double pPhi = tempPos.Lat;
        double pLamb = tempPos.Lon;

        if (Geo2Lamb) {
            Theta = Ntemp * LambDiff(pLamb, Lamb0);

            Ttemp = tnTfunc(pPhi, Ecc);
            Rho = AScaled * Ftemp * pow(Ttemp, Ntemp);

            XYcoord.pXcoord = Rho * Math.sin(Theta);
            XYcoord.pYcoord = Rho0 - Rho * Math.cos(Theta);
        } else {
            if (Ntemp < 0.0) {
                SignN = -1.0;
            } else {
                SignN = 1.0;
            }

            Theta = Math.atan((SignN * XYcoord.pXcoord) / ((SignN * Rho0) - (SignN * XYcoord.pYcoord)));

            TempDouble = Rho0 - XYcoord.pYcoord;
            Rho = Math.sqrt(XYcoord.pXcoord * XYcoord.pXcoord + TempDouble * TempDouble);

            // Rho must take sign of Ntemp
            if ((Rho < 0.0 && Ntemp >= 0.0) || (Rho >= 0.0 && Ntemp < 0.0))
                Rho = -Rho;

            Ttemp = pow(Rho / (AScaled * Ftemp), 1.0 / Ntemp);
            Chi = HALFPI - 2.0 * Math.atan(Ttemp);
            tempPos.Lat = Chi2Phi(Chi, pCurDatum);
            tempPos.Lon = Theta / Ntemp + Lamb0;
        }
    }

    //Stereograph
    boolean UserLL2XYGeo2Ster(double lat, double lon, double lat_c, double k_0,
                              mgnDGeodeticPos3D_t tempPos, XYcoord XYcoord,
                              double false_e, double false_n, double unitfactor,
                              mgnDatumInfo_t localDatum) {

        LoalToDatum(true, localDatum, tempPos);

        Geo2Ster(true, lat, lon, lat_c, k_0, tempPos, XYcoord, localDatum);
        XYcoord.pXcoord = (XYcoord.pXcoord + false_e) / unitfactor;
        XYcoord.pYcoord = (XYcoord.pYcoord + false_n) / unitfactor;
        if ((XYcoord.pXcoord < -9999999.0) || (XYcoord.pXcoord > 99999999.0) ||
                (XYcoord.pYcoord < -9999999.0) || (XYcoord.pYcoord > 99999999.0)) {
            return false;
        }

        return true;
    }

    void Geo2Ster(boolean Geo2Ster,
                  double lat, double lon, double lat_c, double k_0,
                  mgnDGeodeticPos3D_t tempPos, XYcoord XYcoord,
                  mgnDatumInfo_t localDatum) {

        double mgnphi_0 = lat;
        double mgnlamb_0 = lon;
        double mgnphi_c = lat_c;
        double mgnk_0 = k_0;

        mgnGeo2Ster(Geo2Ster, mgnphi_0, mgnlamb_0, mgnphi_c, mgnk_0, tempPos, XYcoord, localDatum);
    }

    double mgnChiFunc(double phi, double ecc) {
        double sin_phi, e_sin_phi;
        // conformal lat

        sin_phi = Math.sin(phi);
        e_sin_phi = ecc * sin_phi;
        return (2.0 * Math.atan(Math.sqrt((1.0 + sin_phi) / (1.0 - sin_phi)
                * pow((1.0 - e_sin_phi) / (1.0 + e_sin_phi), ecc))) - Math.PI / 2.0);
    }

    void mgnGeo2Ster(boolean geo2Ster, double phi_0, double lamb_0,
                     double phi_c, double k_0, mgnDGeodeticPos3D_t tempPos,
                     XYcoord XYcoord, mgnDatumInfo_t pCurDatum) {
        double phi = tempPos.Lat;
        double lamb = tempPos.Lon;

        double temp_double;
        double one_minus_e2, ecc;
        double e_2, e_4, e_6;
        double e_prime_2, e_factor;
        double a_scaled;
        double t_temp, rho, sign_n, chi;
        double half_pi;
        double a_const, b_const, c_const, d_const;
        double chi_0, cos_chi, cos_chi_0;
        double cos_lml0, sin_chi, sin_chi_0;
        double c_e, sin_c_e, cos_c_e;
        double lamb_minus_lamb0, a_temp;

        // get parameters for currently set ellipsoid and calculate several useful things
        e_2 = 1.0 - (1.0 - pCurDatum.flat) * (1.0 - pCurDatum.flat);
        ecc = Math.sqrt(e_2);
        one_minus_e2 = 1.0 - e_2;
        e_4 = e_2 * e_2;
        e_6 = e_4 * e_2;
        e_prime_2 = e_2 / one_minus_e2;
        e_factor = Math.sqrt(pow(1.0 + ecc, 1.0 + ecc) * pow(1.0 - ecc, 1.0 - ecc));
        a_const = -e_2 * 0.25;
        b_const = -e_4 * 0.046875;
        c_const = -e_6 * 0.01953125;
        d_const = 1.0 + a_const + b_const + c_const;
        half_pi = Math.PI / 2.0;

        a_scaled = pCurDatum.smaxis * k_0;

        if (geo2Ster) {
            if (Math.abs(phi_0) < 1.569) {
                // non-polar aspect
                chi = mgnChiFunc(phi, ecc);
                chi_0 = mgnChiFunc(phi_0, ecc);
                sin_chi = Math.sin(chi);
                cos_chi = Math.cos(chi);
                sin_chi_0 = Math.sin(chi_0);
                cos_chi_0 = Math.cos(chi_0);
                lamb_minus_lamb0 = LambDiff(lamb, lamb_0);
                cos_lml0 = Math.cos(lamb_minus_lamb0);
                a_temp = 2.0 * a_scaled * tnMfunc(phi_0, e_2)
                        / (cos_chi_0 * (1.0 + sin_chi_0 * sin_chi
                        + cos_chi_0 * cos_chi * cos_lml0));
                XYcoord.pXcoord = a_temp * cos_chi * Math.sin(lamb_minus_lamb0);

                if (Math.abs(phi_0) > 0.000001) {
                    // oblique aspect
                    XYcoord.pYcoord = a_temp * (cos_chi_0 * sin_chi - sin_chi_0 * cos_chi * cos_lml0);
                } else {
                    // equatorial aspect
                    a_temp = 2.0 * a_scaled / (1.0 + cos_chi * cos_lml0);
                    XYcoord.pYcoord = a_temp * sin_chi;
                }
            } else {
                // polar aspect
                if (phi_0 > 0.0) {
                    sign_n = 1.0;
                    //*************    test   *****************************
                    //  add 180 degrees to lamb_0
                    //**************************************************

                    lamb_0 = InitPIBounds(lamb_0 + Math.PI);
                } else {
                    // for south polar aspect reverse signs of XYcoord.pXcoord, XYcoord.pYcoord, phi_c, phi, lamb, lamb_0
                    sign_n = -1.0;
                }

                lamb_minus_lamb0 = sign_n * LambDiff(lamb, lamb_0);

                if (Math.abs(phi_c) > 1.569) {
                    // phi_c equals 90 degrees, scale applies to pole
                    rho = 2.0 * a_scaled * tnTfunc((sign_n * phi), ecc) / e_factor;
                } else {
                    // scale specified at latitiude phi_c
                    rho = pCurDatum.smaxis * tnMfunc((sign_n * phi_c), e_2) * tnTfunc((sign_n * phi), ecc)
                            / tnTfunc((sign_n * phi_c), ecc);
                }

                XYcoord.pXcoord = sign_n * (rho * Math.sin(lamb_minus_lamb0));
                XYcoord.pYcoord = sign_n * (-rho * Math.cos(lamb_minus_lamb0));
            }
        } else {
            // stereographic to geodetic
            rho = Math.sqrt((XYcoord.pXcoord * XYcoord.pXcoord) + (XYcoord.pYcoord * XYcoord.pYcoord));

            if (Math.abs(phi_0) < 1.569) {
                // non-polar aspect (oblique and equitorial aspects)
                chi_0 = mgnChiFunc(phi_0, ecc);
                sin_chi_0 = Math.sin(chi_0);
                cos_chi_0 = Math.cos(chi_0);

                c_e = 2.0 * Math.atan(rho * cos_chi_0 / (2.0 * a_scaled * tnMfunc(phi_0, e_2)));
                sin_c_e = Math.sin(c_e);
                cos_c_e = Math.cos(c_e);

                chi = Math.asin(cos_c_e * sin_chi_0 + (XYcoord.pYcoord * sin_c_e * cos_chi_0 / rho));

                phi = Chi2Phi(chi, pCurDatum);

                lamb = lamb_0 + Math.atan(XYcoord.pXcoord * sin_c_e /
                        (rho * cos_chi_0 * cos_c_e - XYcoord.pYcoord * sin_chi_0 * sin_c_e));
            } else {
                // polar aspects
                if (phi_0 > 0.0) {
                    sign_n = 1.0;

                    //*************    test   *****************************
                    //  add 180 degrees to lamb_0
                    //**************************************************

                    lamb_0 = InitPIBounds(lamb_0 + Math.PI);
                } else {
                    // for south polar aspect reverse signs of XYcoord.pXcoord, XYcoord.pYcoord, phi_c, phi, lamb, lamb_0
                    sign_n = -1.0;
                }

                if (Math.abs(phi_c) > 1.569) {
                    // phi_c equals 90 degrees, scale applies to pole
                    t_temp = rho * e_factor / (2.0 * a_scaled);
                } else {
                    t_temp = rho * tnTfunc((sign_n * phi_c), ecc)
                            / (pCurDatum.smaxis * tnMfunc((sign_n * phi_c), e_2));
                }

                chi = half_pi - 2.0 * Math.atan(t_temp);

                phi = sign_n * Chi2Phi(chi, pCurDatum);

                temp_double = sign_n * (sign_n * lamb_0 + Math.atan2(sign_n * XYcoord.pXcoord,
                        (-sign_n * XYcoord.pYcoord)));

                // keep -pi < lon <= +pi

                if (temp_double < -Math.PI) {
                    lamb = temp_double + TWOPI;
                } else if (temp_double > Math.PI) {
                    lamb = temp_double - TWOPI;
                } else {
                    lamb = temp_double;
                }
            }
        }

        tempPos.Lat = phi;
        tempPos.Lon = lamb;
    }

    //Polyconic
    boolean UserLL2XYGeo2Poly(double lat, double lon, double k_0,
                              mgnDGeodeticPos3D_t tempPos, XYcoord XYcoord,
                              double false_e, double false_n, double unitfactor,
                              mgnDatumInfo_t localDatum) {

        LoalToDatum(true, localDatum, tempPos);
        Geo2Poly(true, lat, lon, k_0, tempPos, XYcoord, localDatum);

        XYcoord.pXcoord = (XYcoord.pXcoord + false_e) / unitfactor;
        XYcoord.pYcoord = (XYcoord.pYcoord + false_n) / unitfactor;
        if ((XYcoord.pXcoord < -9999999.0) || (XYcoord.pXcoord > 99999999.0) ||
                (XYcoord.pYcoord < -9999999.0) || (XYcoord.pYcoord > 99999999.0)) {
            return false;
        }
        return true;
    }

    void Geo2Poly(boolean Geo2Poly,
                  double lat, double lon, double k_0,
                  mgnDGeodeticPos3D_t tempPos, XYcoord XYcoord, mgnDatumInfo_t localDatum) {
        double mgnlat = lat;
        double mgnlon = lon;
        double mgnk_0 = k_0;
        mgnGeo2Poly(Geo2Poly, mgnlat, mgnlon, mgnk_0, tempPos, XYcoord, localDatum);
    }

    void mgnGeo2Poly(boolean geo2Poly, double phi_0, double lamb_0, double k_0,
                     mgnDGeodeticPos3D_t tempPos,
                     XYcoord XYcoord, mgnDatumInfo_t pCurDatum) {
        double phi = tempPos.Lat;
        double lamb = tempPos.Lon;

        double temp_double;
        double one_minus_e2;
        double e_2, e_4, e_6;
        double e_prime_2;
        double n_temp, m_0;
        double a_scaled;
        double sin_phi, cos_phi, a_scaled_inv, phi_n, c_temp, m_n, cos_2_phi, sin_2_phi;
        double lamb_minus_lamb0, a_temp, b_temp, e_temp, m_n_prime, m_a, m_a_2, phi_n_inc;
        int counter;
        mgnDistArg_t m_args = new mgnDistArg_t();

        // get parameters for currently set ellipsoid and calculate several useful things

        e_2 = 1.0 - (1.0 - pCurDatum.flat) * (1.0 - pCurDatum.flat);
        one_minus_e2 = 1.0 - e_2;
        e_4 = e_2 * e_2;
        e_6 = e_4 * e_2;
        e_prime_2 = e_2 / one_minus_e2;
        m_args.a_const = -e_2 * 0.25;
        m_args.b_const = -e_4 * 0.046875;
        m_args.c_const = -e_6 * 0.01953125;
        m_args.d_const = 1.0 + m_args.a_const + m_args.b_const + m_args.c_const;
        m_args.smaxis = pCurDatum.smaxis;

        //************************ check use of scale factor ****************

        a_scaled = pCurDatum.smaxis * k_0;

        m_args.phi = phi_0;
        m_0 = tnMDist(m_args);

        if (geo2Poly) {
            lamb_minus_lamb0 = LambDiff(lamb, lamb_0);
            if (Math.abs(phi) < 0.000001) {
                XYcoord.pXcoord = a_scaled * lamb_minus_lamb0;
                XYcoord.pYcoord = -m_0 * k_0;
            } else {
                sin_phi = Math.sin(phi);
                cos_phi = Math.cos(phi);
                n_temp = (a_scaled * cos_phi) / (Math.sqrt(1.0 - e_2 * sin_phi * sin_phi) * sin_phi);
                e_temp = lamb_minus_lamb0 * sin_phi;
                XYcoord.pXcoord = n_temp * Math.sin(e_temp);
                m_args.phi = phi;
                XYcoord.pYcoord = tnMDist(m_args) - m_0 + n_temp * (1.0 - Math.cos(e_temp));
            }
        } else {
            a_scaled_inv = 1.0 / a_scaled;
            if ((XYcoord.pYcoord + m_0) < 0.0000001) {
                phi = 0.0;
                lamb = (XYcoord.pXcoord * a_scaled_inv) + lamb_0;
            } else {
                a_temp = (m_0 + XYcoord.pYcoord) * a_scaled_inv;
                b_temp = ((XYcoord.pXcoord * XYcoord.pXcoord) * (a_scaled_inv * a_scaled_inv))
                        + (a_temp * a_temp);
                phi_n = a_temp;

                counter = 0;
                do {
                    ++counter;
                    sin_phi = Math.sin(phi_n);
                    c_temp = Math.sqrt(1.0 - e_2 * sin_phi * sin_phi) * Math.tan(phi_n);
                    m_args.phi = phi_n;
                    m_n = tnMDist(m_args);
                    cos_2_phi = Math.cos(2.0 * phi_n);
                    sin_2_phi = Math.sin(2.0 * phi_n);
                    m_args.a_const = e_2 * 0.25;
                    m_args.b_const = e_4 * 0.09375;
                    m_args.c_const = e_6 * 0.0390625;
                    m_n_prime = -(m_args.a_const + 3.0 * m_args.b_const + 5.0 * m_args.c_const - 1.0) +
                            cos_2_phi * (-3.0 * m_args.a_const - 2.0 * m_args.b_const +
                                    3.0 * m_args.c_const +
                                    cos_2_phi * (5.0 * m_args.b_const + 9.0 * m_args.c_const +
                                            cos_2_phi * (7.0 * m_args.c_const)));
                    m_a = m_n * a_scaled_inv;
                    m_a_2 = m_a * m_a;
                    phi_n_inc = phi_n - (a_temp * (c_temp * m_a + 1.0) - m_a - 0.5 * (m_a_2 + b_temp) * c_temp)
                            / (e_2 * sin_2_phi * (m_a_2 + b_temp - 2.0 * a_temp * m_a) / (4.0 * c_temp)
                            + (a_temp - m_a) * (c_temp * m_n_prime - 2.0 / sin_2_phi) - m_n_prime);

                    temp_double = phi_n;
                    phi_n = phi_n_inc;

                    // converge to about 1/200 second (approx. 0.15 meter).
                    // if longitude of position is greater than 90 degrees from lamb_0,
                    // will not converge.
                } while ((Math.abs(phi_n_inc - temp_double) > 2.424e-8) && (counter < 50));

                phi = phi_n_inc;
                lamb = (Math.asin(XYcoord.pXcoord * c_temp * a_scaled_inv)) / sin_phi + lamb_0;
            }
        }

        tempPos.Lat = phi;
        tempPos.Lon = lamb;
    }

    //Oblique Mercator
    boolean UserLL2XYGeo2Obmer(int two_point, double lat_0,
                               double lon_c, double alpha_c,
                               double lat1, double lon1,
                               double lat2, double lon2, double k_0,
                               mgnDGeodeticPos3D_t tempPos, XYcoord XYcoord,
                               double false_e, double false_n, double unitfactor,
                               mgnDatumInfo_t localDatum) {
        LoalToDatum(true, localDatum, tempPos);
        Geo2Obmer(true, two_point, lat_0, lon_c, alpha_c, lat1, lon1, lat2, lon2, k_0, tempPos, XYcoord, localDatum);

        XYcoord.pXcoord = (XYcoord.pXcoord + false_e) / unitfactor;
        XYcoord.pYcoord = (XYcoord.pYcoord + false_n) / unitfactor;
        if ((XYcoord.pXcoord < -9999999.0) || (XYcoord.pXcoord > 99999999.0) ||
                (XYcoord.pYcoord < -9999999.0) || (XYcoord.pYcoord > 99999999.0)) {
            return false;
        }

        return true;
    }

    void Geo2Obmer(boolean Geo2Obmer,
                   int two_point, double lat_0,
                   double lon_c, double alpha_c,
                   double lat1, double lon1,
                   double lat2, double lon2, double k_0,
                   mgnDGeodeticPos3D_t tempPos, XYcoord XYcoord, mgnDatumInfo_t localDatum) {
        int mgntwo_point = two_point;
        double mgnphi_0 = lat_0;
        double mgnlamb_c = lon_c;
        double mgnalpha_c = alpha_c;
        double mgnphi_1 = lat1;
        double mgnlamb_1 = lon1;
        double mgnphi_2 = lat2;
        double mgnlamb_2 = lon2;
        double mgnk_0 = k_0;

        mgnGeo2Obmer(Geo2Obmer, mgntwo_point, mgnphi_0, mgnlamb_c, mgnalpha_c, mgnphi_1, mgnlamb_1, mgnphi_2, mgnlamb_2, mgnk_0, tempPos, XYcoord, localDatum);

    }

    void mgnGeo2Obmer(boolean geo2Obmer, int two_point, double phi_0,
                      double lamb_c, double alpha_c,
                      double phi_1, double lamb_1, double phi_2, double lamb_2,
                      double k_0, mgnDGeodeticPos3D_t tempPos,
                      XYcoord XYcoord, mgnDatumInfo_t pCurDatum) {
        double phi = tempPos.Lat;
        double lamb = tempPos.Lon;

        double temp_double;
        double one_minus_e2, ecc;
        double e_2, e_4, e_6;
        double e_prime_2;
        double temp_double_1;
        double a_scaled, t_0, t_1, t_2, f_temp;
        double t_temp, sign_n, chi;
        double half_pi;
        double lamb_0;
        double a_const, b_const, c_const, d_const;
        double cos_phi_0, sin_phi_0;
        double a_temp, b_temp, d_temp, e_temp, h_temp, l_temp, g_temp, j_temp, p_temp;
        double l1_minus_l0, gamma_0, sin_gamma_0, cos_gamma_0;
        double v_coord, u_coord, q_temp, s_temp, l_minus_l0;
        double v_temp, u_temp, cos_alpha_c, sin_alpha_c, big_t;

        // get parameters for currently set ellipsoid and calculate several useful things

        e_2 = 1.0 - (1.0 - pCurDatum.flat) * (1.0 - pCurDatum.flat);
        ecc = Math.sqrt(e_2);
        one_minus_e2 = 1.0 - e_2;
        e_4 = e_2 * e_2;
        e_6 = e_4 * e_2;
        e_prime_2 = e_2 / one_minus_e2;
        a_const = -e_2 * 0.25;
        b_const = -e_4 * 0.046875;
        c_const = -e_6 * 0.01953125;
        d_const = 1.0 + a_const + b_const + c_const;
        half_pi = Math.PI / 2.0;

/*******************forward point/az test  *******************************

 geo_2_obmer = 1;
 two_point = 0;
 k_0 = 1.0;
 phi_0 = 36.0 * DEGTORAD;
 lamb_c = -77.7610558 * DEGTORAD;
 alpha_c = 14.3394883 * DEGTORAD;
 phi = 38.8092128 * DEGTORAD;
 lamb = -76.8707953 * DEGTORAD;

 ***********************************************************/

/*************  inverse, two-point test  *******************************

 geo_2_obmer = 0;
 two_point = 1;
 k_0 = 0.9996;
 phi_0 = 40.0 * DEGTORAD;
 phi_1 = 47.5 * DEGTORAD;
 lamb_1 = -122.3 * DEGTORAD;
 phi_2 = 25.7 * DEGTORAD;
 lamb_2 = -80.2 * DEGTORAD;
 XYcoord.pXcoord = 963436.1;
 XYcoord.pYcoord = 4369142.8;

 ***********************************************************/

        a_scaled = pCurDatum.smaxis * k_0;

        cos_phi_0 = Math.cos(phi_0);
        temp_double = cos_phi_0 * cos_phi_0;
        b_temp = Math.sqrt(1.0 + e_2 * temp_double * temp_double / one_minus_e2);
        sin_phi_0 = Math.sin(phi_0);
        temp_double = 1.0 - e_2 * sin_phi_0 * sin_phi_0;
        temp_double_1 = b_temp * Math.sqrt(one_minus_e2);
        a_temp = a_scaled * temp_double_1 / temp_double;
        d_temp = temp_double_1 / (cos_phi_0 * Math.sqrt(temp_double));
        if ((d_temp * d_temp) < 1.0)
            d_temp = 1.0;
        t_0 = tnTfunc(phi_0, ecc);

        if (phi_0 >= 0.0) {
            sign_n = 1.0;
        } else {
            sign_n = -1.0;
        }

        if (two_point > 0) {
            e_temp = (d_temp + sign_n * Math.sqrt((d_temp * d_temp) - 1.0)) * pow(t_0, b_temp);
            t_1 = tnTfunc(phi_1, ecc);
            t_2 = tnTfunc(phi_2, ecc);
            h_temp = pow(t_1, b_temp);
            l_temp = pow(t_2, b_temp);
            f_temp = e_temp / h_temp;
            g_temp = 0.5 * (f_temp - 1.0 / f_temp);
            temp_double = h_temp * l_temp;
            temp_double_1 = e_temp * e_temp;
            j_temp = (temp_double_1 - temp_double) / (temp_double_1 + temp_double);
            p_temp = (l_temp - h_temp) / (l_temp + h_temp);
            temp_double = lamb_1 - lamb_2;
            if (temp_double > Math.PI) {
                lamb_2 += TWOPI;
            } else if (temp_double < -Math.PI) {
                lamb_2 -= TWOPI;
            }
            lamb_0 = LambDiff((0.5 * (lamb_1 + lamb_2)),
                    (Math.atan(j_temp * Math.tan(b_temp * LambDiff(lamb_1, lamb_2) * 0.5) / p_temp) / b_temp));
            l1_minus_l0 = LambDiff(lamb_1, lamb_0);
            gamma_0 = Math.atan(Math.sin(b_temp * l1_minus_l0) / g_temp);
            sin_gamma_0 = Math.sin(gamma_0);
            alpha_c = Math.asin(d_temp * sin_gamma_0);
        } else {
            // point and azimuth define central line
            f_temp = d_temp + sign_n * Math.sqrt(d_temp * d_temp - 1.0);
            e_temp = f_temp * pow(t_0, b_temp);
            g_temp = 0.5 * (f_temp - 1.0 / f_temp);
            gamma_0 = Math.asin(Math.sin(alpha_c) / d_temp);
            sin_gamma_0 = Math.sin(gamma_0);
            lamb_0 = LambDiff(lamb_c, (Math.asin(g_temp * Math.tan(gamma_0)) / b_temp));
        }
        cos_gamma_0 = Math.cos(gamma_0);
        cos_alpha_c = Math.cos(alpha_c);
        sin_alpha_c = Math.sin(alpha_c);

        if (geo2Obmer) {
            if (Math.abs(phi) > (half_pi - 0.0000001)) {
                // phi = +/- pi/2
                if (phi > 0.0) {
                    sign_n = -1.0;
                } else {
                    sign_n = 1.0;
                }
                v_coord = (a_temp / b_temp) * Math.log(Math.tan(0.25 * Math.PI + sign_n * 0.5 * gamma_0));
                u_coord = a_temp * phi / b_temp;
            } else {
                t_temp = tnTfunc(phi, ecc);
                q_temp = e_temp / (pow(t_temp, b_temp));
                s_temp = (q_temp - 1.0 / q_temp) / 2.0;
                big_t = (q_temp + 1.0 / q_temp) / 2.0;
                l_minus_l0 = LambDiff(lamb, lamb_0);
                v_temp = Math.sin(b_temp * l_minus_l0);
                u_temp = (-v_temp * cos_gamma_0 + s_temp * sin_gamma_0) / big_t;

                // if u_temp = +/- 1 then v_coord is infinite.  set v_coord to arbitrary large value

                if (Math.abs(1.0 - Math.abs(u_temp)) < 0.0000001) {
                    v_coord = 99999999.0;
                } else {
                    v_coord = a_temp * Math.log((1.0 - u_temp) / (1.0 + u_temp)) / (2.0 * b_temp);
                }
                temp_double = Math.cos(b_temp * l_minus_l0);
                if (Math.abs(temp_double) < 0.0000001) {
                    u_coord = a_temp * b_temp * l_minus_l0;
                } else {
                    u_coord = a_temp * Math.atan((s_temp * cos_gamma_0 + v_temp * sin_gamma_0)
                            / temp_double) / b_temp;
                }
            }
            XYcoord.pXcoord = v_coord * cos_alpha_c + u_coord * sin_alpha_c;
            XYcoord.pYcoord = u_coord * cos_alpha_c - v_coord * sin_alpha_c;
        } else {
            // oblique mercator to geodetic
            v_coord = XYcoord.pXcoord * cos_alpha_c - XYcoord.pYcoord * sin_alpha_c;
            u_coord = XYcoord.pYcoord * cos_alpha_c + XYcoord.pXcoord * sin_alpha_c;

            q_temp = Math.exp(-b_temp * v_coord / a_temp);
            s_temp = 0.5 * (q_temp - 1.0 / q_temp);
            v_temp = Math.sin(b_temp * u_coord / a_temp);
            u_temp = (v_temp * cos_gamma_0 + s_temp * sin_gamma_0) / (0.5 * (q_temp + 1.0 / q_temp));
            if (Math.abs(1.0 - Math.abs(u_temp)) < 0.0000001) {
                if (u_temp > 0.0) {
                    sign_n = 1.0;
                } else {
                    sign_n = -1.0;
                }
                phi = sign_n * half_pi;
                lamb = lamb_0;
            } else {
                t_temp = pow(e_temp / Math.sqrt((1.0 + u_temp) / (1.0 - u_temp)), 1.0 / b_temp);
                double dTemp = Math.atan(t_temp);
// 			if (dTemp < 0)
// 			{
// 				dTemp += PI;
// 			}
                chi = half_pi - 2.0 * dTemp;
                phi = Chi2Phi(chi, pCurDatum);
                lamb = LambDiff(lamb_0, Math.atan2(s_temp * cos_gamma_0 - v_temp * sin_gamma_0,
                        Math.cos(b_temp * u_coord / a_temp)) / b_temp);
            }
        }

        tempPos.Lat = phi;
        tempPos.Lon = lamb;
    }

    //Britn Grid
    private final static double BRIT_PHI_0 = 49.0 * (Math.PI / 180.0);
    private final static double BRIT_LAMB_0 = -2.0 * (Math.PI / 180.0);
    private final static double BRIT_K_0 = 0.9996012717;
    private final static long BRIT_FALSE_E = 400000L;
    private final static long BRIT_FALSE_N = -100000L;
    private final static double BRIT_E_MIN = 0;
    private final static double BRIT_E_MAX = 1000000;
    private final static double BRIT_N_MIN = 0;
    private final static double BRIT_N_MAX = 1500000;

    public static final String BritLetter[][] = new String[][]{
            {"V", "W", "X", "Y", "Z"},
            {"Q", "R", "S", "T", "U"},
            {"L", "M", "N", "O", "P"},
            {"F", "G", "H", "J", "K"},
            {"A", "B", "C", "D", "E"}
    };

    private boolean GetBritCoord(mgnDGeodeticPos3D_t tempPos, mgnDatumInfo_t localDatum, UTMXY utmxy) {
        LoalToDatum(true, localDatum, tempPos);

        XYcoord XYcoord = new XYcoord();
        XYcoord.pXcoord = utmxy.pXcoord;
        XYcoord.pYcoord = utmxy.pYcoord;
        Geo2Tm(true, BRIT_PHI_0, BRIT_LAMB_0, BRIT_K_0, tempPos, XYcoord, localDatum);

        XYcoord.pXcoord += BRIT_FALSE_E + 0.2;
        XYcoord.pYcoord += BRIT_FALSE_N + 0.2;

        if ((XYcoord.pXcoord < BRIT_E_MIN) || (XYcoord.pXcoord > BRIT_E_MAX) ||
                (XYcoord.pYcoord < BRIT_N_MIN) || (XYcoord.pYcoord > BRIT_N_MAX)) {
            return false;
        }

        int East500kIndex = (int) (XYcoord.pXcoord * 2.0e-6);
        int North500kIndex = (int) (XYcoord.pYcoord * 2.0e-6);
        utmxy.UTMZoneStr = BritLetter[East500kIndex + 1][North500kIndex + 2];

        XYcoord.pXcoord -= East500kIndex * 5.0e+5;
        XYcoord.pYcoord -= North500kIndex * 5.0e+5;

        East500kIndex = (int) (XYcoord.pXcoord * 1.0e-5);
        North500kIndex = (int) (XYcoord.pYcoord * 1.0e-5);
        utmxy.UTMZoneStr += BritLetter[East500kIndex][North500kIndex];

        XYcoord.pXcoord -= East500kIndex * 1.0e+5;
        XYcoord.pYcoord -= North500kIndex * 1.0e+5;

        utmxy.pXcoord = XYcoord.pXcoord;
        utmxy.pYcoord = XYcoord.pYcoord;

        return true;
    }

    int GetBritIndices(String brit_char, ENIndex enIndex) {
        for (enIndex.NorthIndex = 0; enIndex.NorthIndex < 5; ++enIndex.NorthIndex) {
            for (enIndex.EastIndex = 0; enIndex.EastIndex < 5; ++enIndex.EastIndex) {
                if (brit_char.startsWith(BritLetter[enIndex.NorthIndex][enIndex.EastIndex])) {
                    return 1;
                }
            }
        }
        return 0;
    }

    boolean EntBritCoord(mgnDGeodeticPos3D_t tempPos, mgnDatumInfo_t localDatum, UTMXY utmxy) {
        ENIndex EN500kIndex = new ENIndex();
        ENIndex ENIndex = new ENIndex();
        if (GetBritIndices(utmxy.UTMZoneStr, EN500kIndex) == 0 || (EN500kIndex.EastIndex < 2)
                || (EN500kIndex.EastIndex > 3) || (EN500kIndex.NorthIndex < 1) || (EN500kIndex.NorthIndex > 3)
                || GetBritIndices(utmxy.UTMZoneStr.substring(1), ENIndex) == 0) {
            return false;
        }

        utmxy.pXcoord = (EN500kIndex.EastIndex - 2) * 500000L+ ENIndex.EastIndex * 100000L + utmxy.pXcoord - BRIT_FALSE_E;
        utmxy.pYcoord = (EN500kIndex.NorthIndex - 1) * 500000L+ ENIndex.NorthIndex * 100000L + utmxy.pYcoord - BRIT_FALSE_N;

        XYcoord XYcoord = new XYcoord();
        XYcoord.pXcoord = utmxy.pXcoord;
        XYcoord.pYcoord = utmxy.pYcoord;
        Geo2Tm(false, BRIT_PHI_0, BRIT_LAMB_0, BRIT_K_0, tempPos, XYcoord, localDatum);

        LoalToDatum(false, localDatum, tempPos);
        
        return true;
    }

    //大地坐标 => 经纬度
    boolean UTM2Geo(double Xcoord, double Ycoord, mgnDatumInfo_t localDatum, int Zone,
                    mgnDGeodeticPos3D_t tempPos, String letterZone) {

        mgnEntUtmCoord(Xcoord, Ycoord, localDatum, Zone, tempPos, letterZone);
        LoalToDatum(false, localDatum, tempPos);
        return true;
    }

    boolean mgnEntUtmCoord(double Xcoord, double Ycoord, mgnDatumInfo_t pDatumBlock, int Zone,
                           mgnDGeodeticPos3D_t tempPos, String letterZone) {

        XYcoord XYcoord = new XYcoord();
        XYcoord.pXcoord = Xcoord - UTM_FALSE_E;
        XYcoord.pYcoord = Ycoord;
        if (letterZone == "S") {
            XYcoord.pYcoord = Ycoord - UTM_FALSE_N;
        }
//        int zone = mgnUTMZone(Pos.Lat, Pos.Lon);
        mgnGeo2Tm(false, 0, UTMcnmer(Zone), UTM_K_0, tempPos, XYcoord, pDatumBlock);

        return true;
    }

    boolean UserXY2LLGeo2Tm(double lat, double lon, double k_0,
                            mgnDGeodeticPos3D_t tempPos,
                            XYcoord XYcoord,
                            double false_e,
                            double false_n,
                            double unitfactor,
                            mgnDatumInfo_t pDatumBlock) {
        XYcoord.pXcoord = XYcoord.pXcoord * unitfactor - false_e;
        XYcoord.pYcoord = XYcoord.pYcoord * unitfactor - false_n;
        Geo2Tm(false, lat, lon, k_0, tempPos, XYcoord, pDatumBlock);
        LoalToDatum(false, pDatumBlock, tempPos);
        return true;
    }


    boolean MGRS2Geo(XYcoord XYcoord, int nGridScale, String MGRSZoneStr, int nOption, mgnDatumInfo_t pDatumBlock, int Zone,
                     mgnDGeodeticPos3D_t tempPos) {
        for (int i = 0; i < nGridScale; ++i) {
            XYcoord.pXcoord = XYcoord.pXcoord * 10;
            XYcoord.pYcoord = XYcoord.pYcoord * 10;
        }
        return mgnEntMgrsCoord(XYcoord, MGRSZoneStr, nOption, pDatumBlock, Zone, tempPos);
    }

    public static final String MGRS_NS_ZONE_LETTER[] = new String[]{
            "C", "D", "E", "F", "G",
            "H", "J", "K", "L", "M",
            "N", "P", "Q", "R", "S",
            "T", "U", "V", "W", "X"
    };

    boolean mgnEntMgrsCoord(XYcoord XYcoord, String GridStr, int MgrsOpt,
                            mgnDatumInfo_t pDatumBlock, int MgrsZone, mgnDGeodeticPos3D_t tempPos)

    {
        boolean ret_val;
//	mgnC08_t		temp_grid_str1[16], temp_grid_str2[16];
//	mgnC08_t		check_string[6];
        int table_option;
        int i, bad_letter, southern_hem;
        int east_index, north_index;
        int itemp_lat;
        double x_coord, y_coord;
        double temp_double, zone_northing = 0;
        double temp_double_2, temp_double_3;

	/* Check grid letter    */

	    if (GridStr.length() < 1) {
            GridStr = "       ";//死机更改
        }

        north_index = 0;
        for (i = 0; i < 20; ++i) {
            if (GridStr.startsWith(MGRS_NS_ZONE_LETTER[i])) {
                north_index = i;
                break;
            }
        }
        bad_letter = 1;
        for (i = 0; i < 24; ++i) {
            if (GridStr.substring(1).startsWith(MGRS_GRID_LETTER[i])) {
                bad_letter = 0;
                break;
            }
        }
        if (bad_letter == 0) {
            bad_letter = 1;
            for (i = 0; i < 24; ++i) {
                if (GridStr.substring(2).startsWith(MGRS_GRID_LETTER[i])) {
                    bad_letter = 0;
                    break;/*>>>>>>>>>>>>>>>>*/
                }
            }
        }
        if ((MgrsZone < 1) || (MgrsZone > 60) || (north_index == 0) || (bad_letter == 1)) {
            ret_val = false;
        } else {

	/* Convert to latitude and longitude and store  */

            String es = GridStr.substring(1, 2);
            east_index = -1;
            do {
                ++east_index;
            } while (es.compareToIgnoreCase(MGRS_EAST_LETTER2[east_index]) != 0);


            x_coord = (double) ((east_index - 4 - ((MgrsZone % 3) * 8)) * 100000L);
            x_coord += XYcoord.pXcoord;

            itemp_lat = (north_index - 10) * 8;
            temp_double = (itemp_lat * DEGTORAD);
            temp_double_2 = UTMcnmer(MgrsZone);

            mgnDGeodeticPos3D_t tempPos1 = new mgnDGeodeticPos3D_t();
            tempPos1.Lat = temp_double;
            tempPos1.Lon = temp_double_2;
            XYcoord XYcoord1 = new XYcoord();
            mgnGeo2Tm(true, UTM_PHI_0, UTMcnmer(MgrsZone), UTM_K_0,
                    tempPos1, XYcoord1, pDatumBlock);

            zone_northing = XYcoord1.pYcoord;

            southern_hem = 0;
            if (zone_northing < 0.0) {
                zone_northing += 1.0e+7;
                southern_hem = 1;
            }

            String gs = GridStr.substring(2);
            north_index = -1;
            do {
                ++north_index;
            } while (gs.compareToIgnoreCase(MGRS_GRID_LETTER[north_index]) != 0);

            switch (MgrsOpt) {
                case 0: // MGRS_TP:
                    table_option = 1;
                    break;
                case 1: // MGRS_TE:
                    table_option = 2;
                    break;
                default:
                    table_option = 4;
            }

            if ((pDatumBlock.flags & table_option) > 0) {
                if ((MgrsZone % 2) == 0) north_index -= 5;
            } else {
                if (MgrsZone % 2 != 0) {
                    north_index -= 10;
                } else {
                    north_index -= 15;
                }
            }
            if (north_index < 0) north_index += 20;

            y_coord = (double) (north_index * 100000L);

            while ((y_coord < zone_northing) &&
                    ((zone_northing - y_coord) > 200000.0)) {
                y_coord += 2.0e+6;
            }

            y_coord += XYcoord.pYcoord;
            if (southern_hem > 0) y_coord -= 1.0e+7;

//		x_coord = x_coord + 0.2;
//		y_coord = y_coord + 0.2;

            XYcoord.pYcoord = y_coord;
            XYcoord.pXcoord = x_coord;

            mgnGeo2Tm(false, UTM_PHI_0, UTMcnmer(MgrsZone), UTM_K_0, tempPos,
                    XYcoord, pDatumBlock);

//		fmt_mgrs(temp_lat, temp_lon, temp_grid_str1, temp_grid_str2,
//	              pDatumBlock, mgrsgref, mgrstopt);
//		if (memcmp(temp_grid_str1,check_string,5) != 0) {
//			ret_val = FALSE;
//		} else {
//			*dlat = temp_lat;
//			*dlon = temp_lon;
//		}
            ret_val = true;
        }
        return ret_val;
    }

    boolean UserXY2LLGeo2Lam(boolean OneStandPar,
                             double lat, double lon,
                             double lat1, double lat2, double k_0,
                             mgnDGeodeticPos3D_t tempPos,
                             XYcoord XYcoord,
                             double false_e,
                             double false_n,
                             double unitfactor,
                             mgnDatumInfo_t pDatumBlock) {
        XYcoord.pXcoord = XYcoord.pXcoord * unitfactor - false_e;
        XYcoord.pYcoord = XYcoord.pYcoord * unitfactor - false_n;
        Geo2Lam(false, OneStandPar, lat, lon, lat1, lat2, k_0, tempPos, XYcoord, pDatumBlock);
        LoalToDatum(false, pDatumBlock, tempPos);
        return true;
    }

    boolean UserXY2LLGeo2Ster(double lat, double lon, double lat_c, double k_0,
                              mgnDGeodeticPos3D_t tempPos,
                              XYcoord XYcoord,
                              double false_e,
                              double false_n,
                              double unitfactor,
                              mgnDatumInfo_t pDatumBlock) {
        XYcoord.pXcoord = XYcoord.pXcoord * unitfactor - false_e;
        XYcoord.pYcoord = XYcoord.pYcoord * unitfactor - false_n;

        Geo2Ster(false, lat, lon, lat_c, k_0, tempPos, XYcoord, pDatumBlock);
        LoalToDatum(false, pDatumBlock, tempPos);
        return true;
    }

    boolean UserXY2LLGeo2Poly(double lat, double lon, double k_0,
                              mgnDGeodeticPos3D_t tempPos,
                              XYcoord XYcoord,
                              double false_e,
                              double false_n,
                              double unitfactor,
                              mgnDatumInfo_t pDatumBlock) {
        XYcoord.pXcoord = XYcoord.pXcoord * unitfactor - false_e;
        XYcoord.pYcoord = XYcoord.pYcoord * unitfactor - false_n;
        Geo2Poly(false, lat, lon, k_0, tempPos, XYcoord, pDatumBlock);
        LoalToDatum(false, pDatumBlock, tempPos);
        return true;
    }

    boolean UserXY2LLGeo2Obmer(int two_point, double lat_0,
                               double lon_c, double alpha_c,
                               double lat1, double lon1,
                               double lat2, double lon2, double k_0,
                               mgnDGeodeticPos3D_t tempPos,
                               XYcoord XYcoord,
                               double false_e,
                               double false_n,
                               double unitfactor,
                               mgnDatumInfo_t pDatumBlock) {
        XYcoord.pXcoord = XYcoord.pXcoord * unitfactor - false_e;
        XYcoord.pYcoord = XYcoord.pYcoord * unitfactor - false_n;
        Geo2Obmer(false, two_point, lat_0, lon_c, alpha_c, lat1, lon1, lat2, lon2, k_0, tempPos, XYcoord, pDatumBlock);
        LoalToDatum(false, pDatumBlock, tempPos);
        return true;
    }
}
